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Abstract 

This paper presents a Lie- Trotter splitting for inertial Langevin 
equations (Geometric Langevin Algorithm) and analyzes its long-time 
statistical properties. The splitting is defined as a composition of a 
variational integrator with an Ornstein-Uhlenbeck flow. Assuming the 
exact solution and the splitting are geometrically ergodic, the paper 
proves the discrete invariant measure of the splitting approximates 
the invariant measure of inertial Langevin to within the accuracy of 
the variational integrator in representing the Hamiltonian. In par- 
ticular, if the variational integrator admits no energy error, then the 
method samples the invariant measure of inertial Langevin without 
error. Numerical validation is provided using explicit variational inte- 
grators with first, second, and fourth order accuracy. 
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1 Introduction 



Overview This paper analyzes equilibrium statistical accuracy of discretiza- 
tions of inertial Langevin equations based on variational integrators. Varia- 
tional integrators are time-integrators adapted to the structure of mechanical 
systems [9]. The theory of variational integrators includes discrete analogs 
of the Lagrangian, Noether's theorem, the Euler-Lagrange equations, and 
the Legendre transform. Variational integrators can incorporate holonomic 
constraints (via, e.g., Lagrange multipliers) (26] and multiple time steps to 
obtain so-called asynchronous variational integrators [8]. 

The generalization of variational integrators the paper analyzes are de- 
rived from a Lie- Trotter splitting of inertial Langevin equations into Hamil- 
tonian and Ornstein-Uhlenbeck equations. The integrator is then defined 
by selecting a variational integrator to approximate the Hamiltonian flow 
and using the exact Ornstein-Uhlenbeck flow. Such a generalization of vari- 
ational integrators to inertial Langevin equations will be called a Geometric 
Langevin Algorithm (GLA). 

This type of splitting of inertial Langevin equations is natural, but seems 
to have been only recently introduced in the literature (for molecular dynam- 
ics see [3|25], for dissipative particle dynamics see |H1[T9], and for inertial 
particles see [IS])- This paper is geared towards applications in molecular dy- 
namics where inertial Langevin integrators (including the ones cited above) 
have been based on generalizations of the widely used Stormer-Verlet inte- 
grator. The Stormer-Verlet integrator is attractive for molecular dynamics 
because it is an explicit, symmetric, second-order accurate, variational inte- 
grator for Hamilton's equations. In molecular dynamics it was popularized 
by Loup Verlet in 1967. Other popular generalizations of the Stormer-Verlet 
integrator to inertial Langevin equations include Briinger-Brooks-Karplus 
(BBK) [5], van Gunsteren and Berendsen (vGB) [23], and the Langevin- 
Impulse (LI) methods [20J. The LI method is also based on a splitting of 
inertial Langevin equations, but it is different from the splitting considered 
here. To our knowledge there are few results in the literature which quantify 
the long-time statistical accuracy of the Lie- Trotter splitting considered here. 

GLA is not only quasi-symplectic as defined in RL1 and RL2 of [14], 
but also conformally symplectic, i.e., preserves the precise symplectic area 
change associated to the flow of inertial Langevin processes [12] . One way to 
prove this property is by deriving the scheme from a variational principle and 
analyzing its boundary terms as done in the context of stochastic Hamiltonian 
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systems without dissipation in [3J. 

Organization of the Paper In §|2] the main results of the paper are pre- 
sented, states all of the hypotheses used in the paper. These hypotheses 
are invoked in §Hwhere it is proved that GLA is pathwise convergent on finite 
time intervals (Theorem I2.ip . GLA is geometrically ergodic with respect to 
a nearby invariant measure on infinite time intervals (Theorem I2.2[) . and the 
equilibrium statistical accuracy of GLA is governed by the order of accuracy 
of the variational integrator in representing the Hamiltonian (Theorem 12. 3p . 
In §3 numerical validation is provided. In the Appendix we review some 
basic facts on variational integrators for the reader's convenience. 

Limitations In a nutshell the main result of the paper states that if GLA 
is geometrically ergodic with respect to a unique invariant measure, the error 
in sampling the invariant measure of the SDE is determined by the energy 
error in GLA's variational integrator. Now if the inertial Langevin equations 
have nonglobally Lipschitz drift and the GLA is based on an explicit varia- 
tional integrator, GLA may fail to be geometrically ergodic. In particular, 
for any step-size there will be regions in phase space where the Lipschitz 
constant of the drift is beyond the linear stability threshold of GLA's under- 
lying variational integrator. Hence, an explicit GLA will be stochastically 
unstable. Since our results rely on a strong form of stochastic stability of 
GLA (namely, geometric ergodicity), they may not hold in this case. 

To stochastically stabilize GLA, one can use GLA as a proposal move in 
a Metropolis-Hasting method. For a numerical analysis of the Metropolis- 
adjusted scheme, the reader is referred to jl]. A difficulty in Metropolizing 
inertial Langevin is that its solution is not reversible. However, the solution 
composed with a momentum flip is reversible. The role of momentum flips 
in Metropolizing Langevin integrators is qualitatively and computationally 
analyzed in [U0EI7] ■ For a quantitative treatment of the role of momentum 
flips in pathwise accuracy the reader is referred to [I]. 

Extension to manifolds For the sake of clarity, the setting of this paper 
is inertial Langevin equations on a flat space, but we stress GLA and its 
properties generalize to manifolds. We refer to Remark 12.11 and to [2] for 
details. 
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2 Main Results of Paper 

Inertial Langevin The setting of the paper is a dissipative stochastic 
Hamiltonian system (as in [2TJI23]) on M. n , with phase space M 2n , and smooth 
Hamilton H G C°°(]R 2 " , M). In terms of which consider the following inertial 
Langevin equations 



where the following matrices have been introduced: 



" o i 


, c = 


"0 


0" 


-I 





I 



Here W is a standard 2n- dimensional Wiener process, or Brownian motion, 
(3 > is a parameter referred to as the inverse temperature, and 7 > is 
referred to as the friction factor. We will often write the continuous solution 
in component form as Y(£) = (Q(t), P(t)) where Q(t) and P(t) represent 
the instantaneous configuration and momentum of the system, respectively. 
We shall assume the Hamiltonian is separable and quadratic in momentum: 



where M is a symmetric positive definite mass matrix and U is a potential 
energy function. Despite the degenerate diffusion in (JTJ, under certain regu- 
larity conditions on U, the solution to this SDE is geometrically ergodic with 
respect to an invariant probability measure \i with the following density [23J: 



where Z = J R2n exp (—(3H(q,p)) dqdp. The invariant measure /1 is known as 
the Boltzmann-Gibbs measure. 




WH(Y)dt - 1 CVH(Y)dt + v/^FiCdW 
x e R 2n 



(1) 



H(q,p) = -p T M- l p + U(q), 



7r(q,p) = Z- 1 exp(-(3H(q,p)) 



(2) 
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Geometric Langevin Algorithm Let N and h be given, set T = Nh and 
tfc = hk for k = 0,...,N. Observe that the conservative part of ([I]) defines 
Hamilton's equations for the Hamiltonian H: 

dY = JVH(Y)dt 
dQ = M~ x Pdt 

(3) 

dP = -VU{Q)dt v ' 

Let hbe a. fixed step-size. We apply a pth-order accurate variational integra- 
tor, 8h '■ IR 2n — > M. 2n , to approximate the Hamiltonian flow of (j3J) (p > 1). The 
nonconservative part of the inertial Langevin equation defines an Ornstein- 
Uhlenbeck process in momentum governed by the following linear SDE: 



or, 

r-l 



dY = -^CVH(Y)dt + v / 27/F T CdW 



or, 

' dQ = 



dP = -^M- l Pdt + y^dW 



(4) 



Reference [16] aptly refers to flU) as a Gaussian SDE since its stationary 
distribution on M 2n is Gaussian in momentum. 

The following stochastic evolution map ip tk +h,t k '■ ^ 2n - > M. 2n defines the 
stochastic flow of (j4]): 

ipt k +h,t h '■ 

(Q,P) (q,e-~< M ~ lh p+ v^l^V^^^^^^fslj , (5) 

with ?/> SjS (a?) = as and for < r < s < t recall the Chapman-Kolmogorov 
identity ipt,s ° i>s,r{x) = for all x e IR 2ra . For the distribution of the 

solution, the stochastic flow will be denoted simply by iph- To make this map 
explicit, let £ ~ A/"(0, 1) and set 



=/3 -1 (J - exp(-2 7 M" 1 /i)) M 
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and define Ah to be the decomposition matrix arising from the Cholesky 
factorization of S^, i.e., A^A^ = E^. In terms of these, introduce the 
following flow map: 

A--(q,p)^(q,e^ M ' lh p + A h ^). (6) 

In distribution (jSj) is identical to (j3J). 

Given X k G M 2n and h, the Geometric Langevin Algorithm (GLA) is 
defined as the following Lie- Trotter splitting integrator for (pQ): 

X k+ i := 9 h o i) tk +h,t k {X k ) (7) 

for k = 0, N -1 with X = x. 

Remark 2.1. Observe that GLA generalizes to inertial Langevin equations 
on a manifold. This generalization is possible because its symplectic compo- 
nent can be defined as a variational integrator for Hamilton's equations on a 
manifold and its Ornstein-Uhlenbeck component can be defined as the solution 
of an SDE on a vector space. This generalization is motivated by molecu- 
lar systems with holonomic constraints. As mentioned in the introduction, 
variational integrators can incorporate holonomic constraints. In the special 
case that the configuration manifold of GLA is compact (e.g., SO (3) ) and the 
potential energy is smooth, then the assumption on the geometric ergodicity 
of GLA is typically satisfied for sufficiently small time-step. 

Given Z k e R 2n and h, let $ h ■ M 2n -»> K 2n denote the exact time-A flow 
of Hamilton's equations ()3]). The Exact Splitting is defined as 

Z k+1 :=$ h °ipt k+ h,t k (Z k ) (8) 

for k = 0, N - 1 with Z = x. 

Properties of GLA The assumptions that appear in the following theo- 
rems are provided in £j3j 

Let E^j-} denote the expectation conditioned on the initial condition 
being x e R 2n . In terms of this notation, we can quantify the strong conver- 
gence of GLA to solution trajectories of inertial Langevin ([1]). The precise 
statement follows 
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Theorem 2.1 (Pathwise Accuracy). Assume \3.1\ and \3.2 . For any T > 0, 

there exist h c > and C(T) > 0, such that for all h < h c , x G M? n , and 
t G [0,T], GLA satisfies 

(E x {\X [m -Y([t/h\h)\ 2 })^ < C(T)(1 + \x\ 2 ) 1/2 h. (9) 

This result is expected because a Lie- Trotter splitting is first-order for 
deterministic ODEs, and the noise in (OQ) is additive. 

Using this pathwise convergence, it is shown that GLA is geometrically 
ergodic with respect to a discrete invariant measure 

Theorem 2.2 (Geometric Ergodicity). Assume [X71 [3~J% and\EIR Then 
GLA is geometrically ergodic with respect to a discrete invariant measure fih 
and the continuous Lyapunov function (cf. Assumption \3. 3\) . That is, there 
exist h c > 0, A > , and C3 > 0, such that for all h < h c and for all k > 2, 

\E X {f(X k )} - ii h {f)\ < C 3 V(x)e- Xkh , V x E ffi 2n , 

and for all test functions satisfying \f(y)\ < C^V{y) for all y e M. 2n . 

We stress this result is a consequence of strong convergence of GLA and 
the assumptions made on the potential energy and variational integrator. 
These assumptions are sufficient, but not necessary to guarantee this result. 

Using geometric ergodicity we can quantify the equilibrium statistical 
accuracy of GLA. If p represents the global accuracy of GLA's underly- 
ing variational integrator, then is in TV distance 0(h p ) away from the 
Boltzmann-Gibbs measure fi. To be precise, the main result of the paper 
states 

Theorem 2.3 (Long-Run Accuracy). Assume{3J\ E3 and\EM Let fi h 

denote the discrete invariant measure of GLA. Then, there exist C > and 
h c > 0, such that for all h < h c , 

I A* - HhWv < ChP . 

There is a stronger argument in [23] based on the Feynman-Kac formula 
that can extend Theorem 12.31 to 

Hf)-M\ <chP, (10) 

for all test functions / G L^(IR 2n ) that are smooth with polynomial growth 
at infinity. The paper proves Theorem 12.31 with a more direct strategy. An 
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important point is that the proof is transparent since it involves a forward 
error analysis and does not rely on knowing the precise form of u^. The 
proof relies on the existence of fih and the nature of the convergence of GLA 
from a nonequilibrium position. Indeed, a backward error analysis of this 
discretization of the SDE (Tj[|) to characterize this invariant measure would be 
substantially more involved. 

Implications As a consequence of the TV error estimate derived in this 
paper, one can control the order of accuracy of fih by controlling the order of 
accuracy of GLA's underlying variational integrator. This is the distinguish- 
ing feature of GLA. Existing theory would indicate the accuracy of is the 
same order as the weak or strong accuracy of GLA. Theorem 12. II states GLA 
is just first-order accurate on solution trajectories. Hence, existing theory 
would suggest that the equilibrium statistical accuracy of GLA is first-order, 
rather than pth-order accurate (where p is the order of accuracy of GLA's 
underlying variational integrator). 

Existing theory would indicate to obtain a higher-order approximation 
of the invariant measure one would require a higher-order approximant to 
SDE ([T]) which entails approximation of multiple n-dimensional stochastic 
integrals per time-step. It is well-known that such higher-order discretiza- 
tions of SDEs are computationally intensive. In contrast, a step of GLA 
requires evaluation of a single, n-dimensional stochastic integral per time- 
step. According to the main result of this paper, the order of accuracy of the 
variational integrator can be used to tune the TV-distance in Theorem 12.31 
to a desired tolerance. 

3 Preliminaries 

The following assumptions on the potential energy, U : Q — > R, will be used 
in this paper. These hypotheses are the same as those made in §7 of [TU] . 

Assumption 3.1 (Assumptions on Potential Energy). The potential energy 
function U G C°°(R n ,R) satisfies: 

Ul) there exists a real constant A > such that 

|W(<Z ) - W( 9l )| < A \q - q x \, V q , q x E W 1 . 
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U2) there exists a real constant A\ > such that 

U{q) > ^(1 + |g| 2 ), Vgel". 

By standard results in stochastic analysis, condition Ul is sufficient to 
guarantee almost sure existence and pathwise uniqueness of a solution to 
([Q). The condition U2 ensures that e~^ H is integrable over IR 2ra , and hence, 
that the Boltzmann-Gibbs measure is a well-defined probability measure. 
Assuming the solution to (pQ) is geometrically ergodic, we will prove in this 
paper that conditions Ul - U2 together with the following assumptions on 
the variational integrator, 6^ : R 2n — > R 2n , are sufficient (but not necessary) 
to guarantee geometric ergodicty of GLA. 

Assumption 3.2 (Assumptions on Variational Integrator). For any t > 
let d t denote the exact Hamiltonian flow of 03) ■ The variational integrator 
9 h : R 2n — > R 2n satisfies the following. 

VI) 6h is the discrete Hamiltonian map of a hyperregular discrete Lagrangian 
L d : W 1 x R n R (cf. §2$ and 

V2) there exist constants B > and h c > 0, such that for any h < h c , 
\9 h {x) - d h (x)\ < B (l + \x\ 2 ) 1/2 h p+1 , V x e R 2n . 

As discussed in Appendix I, the condition VI implies that 8h is symplec- 
tic, and hence, Lebesgue measure preserving. It will also be an important 
ingredient in proving Theorem 12 .21 on geometric ergodicity of GLA. The con- 
dition V2 states that the integrator is locally (p + l)t/i-order accurate. 

Finally, we make the following structural assumption on ([1]). 

Assumption 3.3 (Existence of a Lyapunov Function). There exists V G 
C°°(]R 2n ,]R) and constants d > such that 

C (l + |a;| 2 ) < V(x) < Ci(l + \x\ 2 ), VV(x)<C 2 (l + \x\), V x G R 2n , 

lima-^oo V(x) = oo, a > and c > 0, such that for all t > 0, 

E x {V{Y(t))} < e- at V{x) + -(1 - e~ at ), V i 6 R 2n . 

a 
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4 Analysis of GLA 



4.1 Pathwise Convergence 

Here GLA is shown to be first-order mean-squared convergent, which is a 
notion of pathwise convergence to solutions of ([T]) |15],[22] . The first-order 
accuracy of GLA on solution trajectories is not surprising because the method 
is derived from a Lie- Trotter splitting of ([I]). It is simply a generalization of 
the well-known fact that Lie- Trotter splittings of deterministic ODEs yield 
first-order accurate methods. This generalization is possible despite the lack 
of regularity in solutions because the noise in ([!]) is additive. Since the proof 
is standard, it will be kept terse. 

Theorem 2.1 (Pathwise Accuracy). Assume \3.1\ and \3.2 . For any T > 0, 

there exist h c > and C(T) > 0, such that for all h < h c , x G M? n , and 

te[o,T], 

(E x {\X [m - Y{[t/h\h)\ 2 }fl 2 < C{T){\ + \x\ 2 fl 2 h. (11) 

Proof. By standard results in stochastic analysis, condition Ul guarantees 
there a.s. exists a pathwise unique solution to ([!]): Y(t) G M 2n for t G [0, T] 
with Y(0) = x. Moreover, one can obtain the following bound on the second 
moment of of the solution: for all T > 0, there exists a C(T) > such that 
for all t G [0,T], 

E x {\Y(t)\ 2 } <C(T)(l + \x\ 2 ). (12) 

We will use this bound to invoke Theorem 1.1 in [15] which enables one 
to deduce global mean-squared error estimates of a discretization from lo- 
cal mean-squared error and local mean deviation. First, we establish this 
estimate for the exact splitting OH]). By using Assumption Ul, it is straight- 
forward to show (see Lemma [7. ip that there exists C > such that: 

\E x {Y(h)- Z 1 }\ < C(l+ \x\ 2 ) 1/2 h 2 (13) 

and 

(E x {\Y(h) - Z^ 2 }) 1 ' 2 <C(1 + \x\ 2 ) 1 ' 2 2 (14) 

Together with ( TT2]) this implies there exist h c > and C(T) > 0, such that 
for all h<h c ,te [0, T] and x G R 2n : 

E-{|Z Lt/fcJ | a }<C(T)(l + |«| 2 ) (15) 
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Hence, by Theorem 1.1 in [15], one can show that for all T > 0, there exist 
h e > and C{T) > 0, such that for all h < h c , t G [0, T] and x G R 2n : 

(E-{|Z Lt/fcJ - Y(\t/h\h)\ 2 })^ < C(T) (1 + \x\ 2 ) 1/2 h. (16) 

Observe that the difference between a single step of GLA (J7|) and the 
exact splitting OH]) can be written as 

x 1 -z 1 = (e h -# h )o></> hfi (x). 

Using Assumption V2 one can show there exists C > such that 

(E X {\X 1 - Zxl 2 }) 172 < C (1 + |:z| 2 ) 1/2 (17) 

and, by Jensen's inequality: 

\E X {X 1 - Zi}| < C (1 + \x\ 2 ) 1/2 h p+1 . (18) 

Together with ( Tl5|) this implies that there exist h c > and C(T) > 0, such 
that for all h < h c , t G [0, T] and a; G M 2n : 

E^lXL^jl^^C^a + W 2 ) (19) 

Using Assumption Ul and Theorem 1.1 of [15], one can also show that for 
all T > 0, there exist h c > and CiT) > 0, such that for all h < h c , t G [0, T] 
and x G M 2n : 

(E"{|X Lt/fcJ - Z^jl 2 }) 1 / 2 < C(T) (1 + l^l 2 ) 172 ^. (20) 

In other words, GLA is 0(h p ) strongly convergent to the exact splitting. One 
can then use the triangle inequality to obtain the estimate in the theorem 
from (1201) and (flB]) . i.e., 

(E x {\X [t/h] -Y([t/h\h)\ 2 })^< 

(E*{\X mi - Z^jl 2 }) 1 ^ + (E*{\Z mi - Y{\t/h\h)\ 2 }fl\ 

V v ' S v ' 

<K(T)(l+\x\2) 1/2 hP <K(T)(l+\x\ 2 ) 1/2 h 

In sum, GLA is first-order strongly convergent to solutions of (pQ). □ 
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4.2 Geometric Ergodicity 

Geometric ergodicity is a strong type of stochastic stability of a Markov 
chain [13]. In this section geometric ergodicity of GLA is established following 
the recipe provided in §7 of jlOJ . In the context of this paper, geometric 
ergodicity means, 

Definition 4.1. A Markov chain X k is said to be geometrically ergodic if 
there exist probability measure p^, p < 1, and M G C°°(M. 2n , such that 

\E x {f(X k )}-p 00 (f)\<M(x)p k , \/xeR 2n , VfceN, (21) 

and for all f G L 2 ^ (R 2n ) satisfying \f(y)\ < M(y) for all y G R 2n . 

Under the hypotheses below, the Lyapunov function from Assumption l3.3l 
is inherited by GLA. 

Theorem 2.2 (Geometric Ergodicity). Assume {3J\ [23 and\EM Then 
GLA is geometrically ergodic with respect to a discrete invariant measure ph 
and the continuous Lyapunov function (cf. Assumption \3. 3\) . That is, there 
exist h c > 0, A > , and C3 > ; such that for all h < h c and for all k > 2, 

\E X {f(X k )} - p h (f)\ < C 3 V(x)e- Xkh , V x G R 2n , 

and for all test functions satisfying \f(y)\ < C^V{y) for all y G M. 2n . 

Proof. This proof is an application of Theorem 2.5 of [10] . To invoke this 
theorem, we will show that GLA inherits the Lyapunov function V : M. 2n — >• M 
of the continuous solution (cf. Assumption 13. 3p and satisfies a minorization 
condition when sampled every other step. 

To prove that GLA inherits the Lyapunov function V : M 2n — > M. we 
use Theorem 7.2 of [10]. This theorem assumes that the Lyapunov function 
of the SDE is essentially quadratic which follows from Assumption I3.3[ and 
that the discretization of the SDE satisfies Condition 7.1 of [10J. Condition 
7.1 (i) is a consequence of a single-step mean-squared error estimate of GLA 
which can be derived from f fl4|) and fflTI) . Condition 7.1 (ii) is satisfied for 
the first and second moments of GLA due to the estimate (fT9]) . Hence all 
of the assumptions of Theorem 7.2 [10] are satisfied, and one can conclude 
that GLA inherits the Lyapunov function V : M 2n — > M up to a constant 
pre-factor. 
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Next, we prove that GLA satisfies a minorization condition when sampled 
every other step. This property follows from Lemma 2.3 of [TO], because 
GLA sampled every other step admits a strictly positive, smooth transition 
probability function. In fact, this transition probability qh : M 2 "xf 2n — > [0, 1] 
can be explicitly characterized, and by inspection it is clear that it is smooth 
as a function of its arguments and strictly positive everywhere. 

To derive this expression, let o/, : I" x K™ 4 [0, 1] denote the transition 
probability of the Ornstein-Uhlenbeck flow iph ©. By a change of variables, 
it's transition density is given explicitly by: 

Oh(P ,Pl) = 



(22) 



where 



S ft = /T 1 (Id- exp(-2 7 M" 1 /i)) M. 

Let D = R 2n x R 2n x R 2n x R n . Since the maps 8 h and i) h enjoy the Markov 
property, the transition probability of the composition 9 h o ip h o 9 h o ip h can 
be expressed as a product of the transition probabilities of its components: 

Qh((q,p),(q,P)) = 

o h (p, Pi)S((q 1 ,p 2 ) - 9 h (q, Px))o h (p 2 , p 3 )S((q, p) - 6 h (q 1 ,p 3 ))dp 1 dp 2 dp 3 dq 1 . 



D 



The zero of the argument of the second Dirac-delta measure (from left) occurs 
at (q 1 ,p 3 ) = 6 l / ^ 1 (q, p). Hence, the above expression simplifies, 

Qh((q,p),(q,P)) = 

o h (p,p 1 )S((q 1 ,p 2 ) - fc (g,Pi))o h (p 2 ,p & )dp 1 dp 2 . 

By condition VI on 6h, the zero of the argument of the remaining Dirac-delta 
measure above is uniquely determined by the discrete Hamiltonian flow of 
the discrete Lagrangian (cf. floT]) in Appendix I). Hence, one obtains: 

<lh ((q,p), (q,p)) = 

I det(D 12 L d (q, q v h))\o h (p, -D x L A {q, q lt h))o h (D 2 L d (q, q ± , h),p 3 ) (23) 
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where (<?i,p 3 ) = 6 l ^" 1 (g,p). Using the hyperregularity assumption on the 
variational integrator VI (cf. ( 15^]) ). ( I2"2"j) . and (I2"3"j) . it is clear that qh is a 
smooth probability transition function that is everywhere strictly positive. 
Hence, by Lemma 2.3 of [10], GLA sampled every other step satisfies a mi- 
norization condition. 

In sum, we have shown that GLA satisfies a minorization condition and 
admits a Lyapunov function. The result follows from invoking Theorem 2.5 
in [101. □ 



4.3 Long- Run Accuracy 

Now we quantify the accuracy of GLA in sampling from the equilibrium 
measure of (Q. For this purpose recall the following definition. 

Definition 4.2 (Invariance of Measure). A Markov chain X k G M 2n is said 
to preserve a probability measure /ioo if for all f G L 2 oo (IR 2ri ) and fc£N, 

E flac E x {f(X k )}=^ OQ (f) (24) 

where /ioo(/) = J R 2 n fd^oo and E lloB E x denotes expectation conditioned on 
the initial distribution being sampled from floo, i.e., 

E IMoc E x {f(X k )} = [ E" s {/(X fc )}/i 0O (£te). 

Given a step-size h, define the deviation GLA makes in preserving the 
Boltzmann-Gibbs measure, /i, as : L 2 (R 2n ) — > R: 

A k h (f) :=E,E x {f(X k )}-^(f). 

Observe that if GLA exactly preserves fi then: 

A*(/) = 0, V/gL^R 2 "). 

The following local error result follows from the Ornstein-Uhlenbeck flow iph 
preserving [i and the variational integrator 6h preserving Lebesgue measure. 

Lemma 4.3. Suppose the potential energy satisfies U2. For a given f G 
L 2 (R 2 «) ; 

K(f) = [ fix) (e-^O-WJ-irW) _ A Kdx) , 
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Proof. The condition U2 ensures that fi is a well-defined probability mea- 
sure. According to the definition of GLA (0), Xi = Oh o ij) h (x). Substitute 
this expression into to obtain: 

Al(f)= [ E x {f(9 h oip h (x))}»(dx)- [ fdii. 

JR 2 " JR 2n 

Since iph preserves \x and Oh is deterministic it follows that, 

K(f) = [ f(e h (x))fi(dx) - f fdfi. 

jR 2n Jr 2 ™ 

Changing variables under the map Oh in the first integral above, and using the 
volume-preserving property of the variational integrator Oh (See Appendix.) 
one obtains the desired expression. □ 

Remark 4.1. As a consequence of Lemma \4-3\ if Oh admits no energy error, 
then GLA preserves //. In particular, the exact splitting (jS]) preserves fi. 

In the situation where GLA is geometrically ergodic, this paragraph quan- 
tifies the equilibrium error of GLA in preserving the BG measure. 

Lemma 4.4. Assume [3J[\3.£\ and \3.3[ Then, there exist C > and h c > 0, 
such that for all h < h c , 

lim \AZ{f)\<ChT, 

and for all f e L 2 ^ (R 2n ) satisfying \f(y)\ < C 3 V{y) for all y e R 2n . 

Proof. Let / G L 2 (R 2n ) such that \f(y)\ < C 3 V(y) for all y E R 2n . The 
term E /J E a: {/(X A r)} can be written as a telescoping sum: 

N 

E,E*{f(X N )} = /!(/) + (V{f(X k )} - E^"{f(X k ^)}) . 

k=l 

By Lemma [4. 3[ one can rewrite/reindex this sum as: 
» N-l 

A H(f)= / y)E-{/(X fc )}fc-«^ 1W) - HW) -l)A*(d«)- (25) 
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Since Oh preserves Lebesgue measure, one can write this deviation as: 

KU) = 

r N-l 

v * / o L 



k=0 

Deviation from Equilibrium _ . , r . ,, , _ 

Energy Error ot Variational Integrator 



(26) 



From ( 126]) it is clear that the equilibrium BG error is due to: 1) how fast GLA 
converges to equilibrium and 2) the local accuracy with which 0^ represents 
the Hamiltonian function H. The equality ( |26|) is the crux of the proof, and 
what follows is an approach to bound A^ (/). 

Since GLA is geometrically ergodic (cf. Theorem 12.21) . one can bound 
A^ (/) from above by 

K(/)| < f Ve-*| C 3 / V{x) e ^mo- h \ x ))-H( x) ) _ j ^ dx y 

Changing variables in the right-hand-side under the map Oh, one can rewrite 
this bound as, 

|Af (/)| < f V e~ Xhk ) C 3 [ V(O h (x)) \ e -Pme h { x ))-H{ x )) _ i\^ dx y 
\k=o J J ^ n 

In the limit as N — > oo, the right-hand-side of the above can be written in 
terms of the formula for the geometric series for e~~ Xh : 

lim \A»(f)\ < — %^ [ V(O h (x)) | c H>(ir(M«0)-*(«0) _ l\^(dx). (27) 



Using the natural bound \e x — 1| < e^' — 1 for all x G IR, one can further 
bound |A*(/)| by: 

lim |A?(/)| < — %^ [ V(O h (x)) ^mo h ( X ))-H( x) \ _ 1} {dx)m (2g) 
Introduce the exact flow dh of Hamilton's equations ([3]) into this bound, 

lim \A N h {f)\ < [ V(O h (x)) - 1) n(dx). 

(29) 
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Set y Q = 9h(x) and y 1 = $h( x )- By the fundamental theorem of calculus, 

H{y l ) - H(y ) = f VH(y + sfa - y )) ■ ( Vl - y )ds. 
Jo 

Using condition Ul and the Cauchy-Schwartz inequality, it follows from the 
above that there exists C > such that 

\H{ Vl ) - H(y )\ < C(l + \ Vl \ + |y |)| tfl - y \. 

Another application of the condition Ul and V2 implies there exists C > 
such that 

\H( yi )-H(y )\<C(l + \x\ 2 )h^\ 

Therefore, 

lim |A?(/)| < — %- h I V(6 h (x)) (eWM-n^ 1 _ i) ^ dx) . ( 30 ) 



Now we show how the the factor V($h(x)) above is handled. 

Since the Lyapunov function is quadratically bounded, the variational in- 
tegrator satisfies V2, and the Hamiltonian vector field is uniformly Lipschitz 
by condition Ul, there exists C > such that 

Km < ^ jfjl + laf ) (e— H«„- _ t ) , {dx) , (31) 

By condition U2 the total energy is quadratically bounded from below. Con- 
sequently one can bound e^ 1311 ^ by e~P D ( 1+ \ x \ -* for some constant D > 0. 
Thus, 

lim \A»(f)\ < 

° [ (1 + laf) (eW+M 2 )** 1 _ i) e -^ 1+ ^dx. 
1 — e Jr 2 ™ ^ ' 

When h < h c = (D/ KY^ P+1 ^ the above integral is finite and one obtains the 
desired error estimate. □ 

A simple application of Theorem 12.31 implies an error estimate for For 
this purpose we introduce the total variation between measures \i and v. 



|/i — v\tv — SU P 
l/l<i 



f(x)(fj,(dx) — v{dx)) 
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Since M(y) > 1 for all y G M 2n , Theorem O applies for all / G L^(M 2n ) 
such that |/(y)| < 1 for all y G M 2n . The TV norm can be written as: 



sup 

l/l<i 



fdfi - E^E x {f(X N )} + E At E a3 {/(X w )} - / fdn 



By the triangle inequality, 

\n - ^|tv < sup |A^(/)| + sup \E,E x {f(X N )} - fi h (f)\ . (32) 

l/l<i l/l<i 

However, under the hypotheses of the theorem, GLA is geometrically ergodic 
with respect to fih and hence, 

lim sup |E,F{/(I ff )} - M\ -> (33) 

N ^°° l/|<i 

and, 

\fJ> - y.h\TV < lim sup |Aft(/)| . (34) 

JV-voo i^-l 

Lemma 14.41 can now be invoked to obtain from (1341 an upper bound for the 
TV distance between /i and jih- This concludes the proof of Theorem 12.31 
which we restate: 

Theorem 2.3 (Long-Run Accuracy). Assume \3~1\ Iff.ffi and Iff. 31 Let fi h 

denote the discrete invariant measure of GLA. Then, there exist C > and 
h c > 0, such that for all h < h c , 

| A« - A^Itv < Ch p . 

In summary, the preceding analysis showed the TV error estimate in The- 
orem 12.31 relies on GLA's variational integrator Oh being volume-preserving 
and pth-order accurate, the Ornstein-Uhlenbeck map iph exactly preserving 
the Boltzmann-Gibbs measure, and GLA being geometrically ergodic. To 
establish the latter, we used the strategy adopted in [TO] which relates path- 
wise convergence of a discretization of an SDE to geometric ergodicity of 
the discretization. This strategy requires the potential force is uniformly 
Lipschitz. 
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5 Validation 



This section tests three different instances of GLA on a variety of simple 
mechanical systems governed by Langevin equations. The purpose of this 
section is to confirm the error estimates provided in the paper. 

Let h be a fixed step size and ~ A/"(0, 1) for k G N. The following 
update scheme is obtained by composing the explicit first-order, symplectic 
Euler method with iph'- 



Pk 

Qk+l 

Pk+1 



-7/1 



Pk 



Pk + 
hp k , 



-2~/h 



6 



(35) 



k+l, 



for k G N. The following integrator is obtained by composing the second- 
order accurate explicit, symmetric, symplectic Stormer-Verlet method with 
tph- 



Pk 

pl/2 
k 

Qk+l 
Pk+1 



Pk + 



-2-yh 



P 



€k, 



Pk ~ 



hdU 
2 dq 



q k + hPf 



1/2 



(36) 



P, 



1/2 



hdU( n \ 
k ~ 2 S^Wfe+lJ) 

for k G N. The following integrator is obtained by composing a fourth- 
order accurate explicit, symmetric, symplectic method due to F. Neri (see, 
e.g., [27]) with ip h : 



Qk, 




e lh p k + 



-2-yh 



P% c ih d Q q 



-- Q l 

Q51 



p 

(Qi 

dihPi+i, 



(37) 



for k G N, and where we have introduced the following constants: 

1 1 - 2 1 / 3 



ci = c 4 



d\ = d 3 



2(2 



2V3) : 
1 



c 2 = c 3 



2 - 2 x /3- 



2(2-2V3)- 
-2 1 / 3 

d 4 = 0. 



2-2V3- 
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The purpose of this fourth-order symplectic integrator is for validation. For 
"optimal" fourth and fifth-order accurate symplectic integrators that mini- 
mize the error in the Hamiltonian, the reader is referred to |11J . 

We will show that despite the fact that fl36|) and fl37|) are only first-order 
pathwise convergent according to Theorem 12.11 they approximate ensem- 
ble averages of /x-integrable functions that satisfy \f(q,p)\ < M(q,p) for all 
(q,p) G M. 2n to within second and fourth-order accuracy, respectively. This 
is consistent with Theorem 12.31 



Linear Oscillator This section follows the analysis of numerical methods 
for linear oscillators governed by Langevin equations developed in [6|fT5] . The 
governing equations for a linear oscillator of unit mass at uniform tempera- 
ture 1//3 are given explicitly by evaluating (JTJ at U(q) = q 2 /2: 

dq = pdt, 

dp = -qdt - 7 pdi + y/2p-^dW. 

The resulting process is Gaussian with stationary distribution given by the 
BG distribution: 

P 00 (q,p) = Z- 1 ex V (-P (4 + ^ 



2 2 



and with 



/i(g 2 ) = lim E{qf} = 1/(3, ) = I™ E{pf} = 1//3, K (qp) = lim E{q tPt } = 0. 

t— >oo t— >oo t—^oo 

The stationary distribution of the geometric Langevin integrators ([35]) - 
(|37|) is also Gaussian with equilibrium distribution of the form: 



where 



a 2 k 
k a 2 



o 2 = lim E{q£}, a 2 = lim E{p 2 n }, « = lim E{q n p n }. 
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This stationary correlation matrix can be explicitly determined. For (1351) its 
entries are given by: 

rr 2 ( 1 + e7fe ) 2 1 , ml 

a « = (2 + 2e^-hi)P = l3 +0{h) 

2 2 + 2e^ fe - h 2 + e 2 ^ 2 1 2 
a ?= (2 + 2e^-^)/3 =^ + ^) 

7 (2 + 2e^ - /i 2 ) (3 [ 1 
Observe that the cumulative error (|35|) makes is of O(h), i.e., 

IS 2 " M<? 2 )l + - M?) 2 )l + |« - fi(qp)\ < 0(h). 
Whereas for (|36|) its entries are given by: 

K = 

and its cumulative error is of 0(h 2 ), i.e., 

K - M<? 2 )l + k 2 - Mp) 2 )I + 1* - /*(gp) | < o(^ 2 ). 

For (13"T|) its entries are given by: 



, 1 (-4-3 x ^2-2x2 2 / 3 U 4 
CT ? = ~Q + ~ 7777, — + °( h " 

1 



(3 144/3 



2 

K = 



and its cumulative error is of 0(h A ), i.e., 

K 2 - M<? 2 )l + K - M 2 )\ + 1* - m?p)I < <^ 4 )- 
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Finally, consider the exact splitting applied to the linear oscillator at 
uniform temperature. Hamilton's equations for a linear oscillator are: 



1 
-1 



(0) 



Po 



with explicit solution given by: 

(*) = 

Thus, the exact splitting update is given by: 



cos(i) 


sin(t) 




Qo 


— sin(i) 


cos(t) 




Po 







cos(/i) 


sin(/i) 




9o 


Pk 




— sin(/i) 


cos(/i) 











sin(/i) 
cos(/i) 



Co- 



In this situation one can show there is no error made in the stationary correla- 
tion matrix. This follows from the fact that the exact solution of Hamilton's 
equations is volume and energy preserving. 



Nonglobally Lipschitz, Nonlinear Oscillator The theory in this paper 
does not apply to this example since the potential force is nonglobally Lips- 
chitz. With a nonglobally Lipschitz potential force, for any h > there will 
exist regions in phase space where the Lipschitz constant of the potential 
force is beyond the linear stability threshold of an explicit variational inte- 
grator 9h- Hence, a GLA based on an explicit variational integrator will be 
stochastically unstable; transient, to be precise. However, for the step-sizes 
and variational integrators employed, and for the duration of the numerical 
experiments, discrete orbits of GLA seem to be confined to a compact re- 
gion of phase space where the variational integrator 9h is linearly stable and 
Monte Carlo estimates are consistent with the error estimates in the paper. 

The governing equations for a cubic oscillator of unit mass at uniform 
temperature 1/(3 are given explicitly by evaluating ([TJ at U(q) = g 4 /4 — g 2 /2: 

f dq = pdt, 

I dp = (q- q 3 )dt - jpdt + ^/2f3- 1 1 dW. 
The resulting potential force is only locally Lipschitz. 
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Time-Step Number of Steps (|35|) (|36|) (1371) 



h 


N 


3.11e-02 


8.03e-03 


1.45e-02 


h/2 


2 N 


1.49e-02 


1.94e-03 


9.80e-04 


h/4 


4 N 


7.42e-03 


4.83e-04 


7.35e-05 



h/8 8 N 3.74e-03 1.29e-04 5.79e-06 



Table 1: The table estimates |/ih(g 2 ) — /u(g 2 )| using empirical time-averages 
with = 40 x 10 9 steps and h = 0.4 with GLA as determined by ( 135]) - 
(|37|) . For subsequent rows the time-steps are halved and the number of steps 
doubled, so that the time-interval of integration is fixed for all experiments. 
The results show that as the time-steps are halved the difference decreases 
linearly for (1331) . nearly quadratically for (ESJ), and nearly quartically for (IBTj) . 
These results are consistent with the error estimates in the paper. 

The estimates shown earlier predict that 

IM<? 2 )-M<? 2 )I< o{hP) 

where p is the order of accuracy of Oh- Hence, one expects near fourth- 
order accuracy for (|37l) . near second-order accuracy for ( 1361) and first-order 
accuracy for (1351) as shown in table [TJ The tests will apply (I35l) -(l37"l) to 
estimate 

f_° q 2 e~ l3U{ - q) dq 
lim E{g 2 } = M (g 2 ) = j^-pu^ 

by empirical averages of the form 

jh,N 

As nicely discussed in [23J, in addition to the discretization error \[i{q 2 ) — 
^h{q 2 ) | one has to cope with the statistical error arising from the time-average 
being finite, i.e., I h,N ~ ^h{q 2 )- The computations were performed with 
7 = 1 and an inverse temperature value of (3 = 2. 

6 Conclusion 

The analysis in this paper represents a first step towards a deeper analysis 
of GLA for molecular systems. In this paper we make assumptions on the 
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Hamiltonian that ensure the solution to inertial Langevin and GLA are ge- 
ometrically ergodic. In particular, we assume the Hamiltonian vector field 
is uniformly Lipschitz and the Hamiltonian is coercive. These hypotheses 
are sufficient to ensure GLA is geometrically ergodic whenever the solution 
process is. In particular, the former hypothesis is important to ensure GLA 
is stochastically stable [T3]. If GLA's underlying variational integrator is 
not globally linearly stable, one can show GLA defines a transient Markov 
chain. Still one can use GLA as proposal step within a Metropolis-Hastings 
algorithm to obtain a stochastically stable Metropolis- Adjusted Geometric 
Langevin Algorithm (MAGLA). A numerical analysis of MAGLA including 
pathwise convergence can be found in [I]. 

A closer inspection of the proof of Theorem 12.31 reveals that the estimate 
relies on the following important ingredients: 

1. GLA is geometrically ergodic with respect to a probability measure 

2. the variational integrator is Lebesgue-measure preserving; 

3. the Ornstein-Uhlenbeck flow preserves /x; and, 

4. the local energy error of the variational integrator is (p + l)t/i-order 
accurate. 

Therefore, we stress that the result holds under more general conditions. The 
main point being: 

If GLA is geometrically ergodic with respect to a unique invariant 
measure, the error in sampling the invariant measure of the SDE 
is determined by the energy error in GLA 's variational integrator. 

7 Appendix 

7.1 Single-Step Error 

Lemma 7.1. Assume \3.1\ and\ 3.2[ Fork small enough, there exists a C > 
such that 

\E x {Y(h) - Zi}| < C(l + \x\ 2 ) 1/2 h 2 (40) 

and 

(E x {\Y(h) - Zxl 2 }) 172 < C (1 + |*| 2 ) 1/2 ti" 2 . (41) 
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Proof. Write Y(£) = (Q(t), P(t)) where Q(t) and P(t) represent the in- 
stantaneous configuration and momentum of the system, respectively. In 
terms of which write the SDE ([I]) as: 

idQ/dt = M l P 

\dP = -VU{Q)dt - 1 M~ 1 Pdt + ^/Q^RdW 

Q(0) = Q and -P(O) = P . It will be useful to write out the solution of 
For this purpose integrate (|42p to obtain: 



Q(h) = Q + hM- 1 P Q + / M" 1 [-Vf/(Q(s)) -jM- 1 P(s)](h-s)ds 

Jo 

+ / (h- s)M- l dW{s) (43) 

Jo 



and 



ph T T 

P[h) = e^ M ~ lh P - hVU(Q ) - (h- S )-_(Q( S )) ■ M- 1 P(s)ds 

Jo 

+ / (I - e~ lM ~ 1{h ~ s) )VU(Q(s))ds + rj (44) 
Jo 



where we have introduced: 



rh 

n = v / 2t^ T / e^ M ~ 1{h - s) dW{s) 
Jo 



Write Z(£) = (Q(t), P(t)) where Q{t) and -P(t) represent the instanta- 
neous configuration and momentum of the exact splitting, respectively. The 
exact splitting after a single step solves 

idQ/dt =M' l P 

(dP/dt = -VU{Q) 1 ' 



where Q(0) = Q and P(0) = e" 7M lfl P + 77. Integrating (SSD yield 



s. 



rh 

Q(h) = Q + hM- l e^' M lh P - / M- l VU(Q(s))(h - s)ds + hM^rf 

Jo 

(46) 
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and 

P(h) = e^ M ' lh P - hVU(Q ) 

- J\h - s) d ^{Q{s)) ■ M- 1 P(s)ds + 77. (47) 

To obtain the mean-squared and mean error estimates we will use the 
following bounds on the second moment of the continuous solution and the 
exact splitting. Namely, for all t G [0, h], there exists a C > such that 

E x {\Z(t)\ 2 } VE X {\Y(t)\ 2 } <C(l + |a;| 2 ) (48) 

where x = (Q ,P ). We will prove this estimate for the exact splitting, 
and omit the proof for the continuous solution since it is very similar. Let 
x = (Q(0),P(0)). By Taylor's formula, 

\Z(t)\ 2 = \x\ 2 + 2^ (^Q{s),M' 1 P{s)^ds + 2 J (P{s),-VU{Q{s))^ds 

By Young's inequality, 
\Z(t)\ 2 < \x\ 2 

+ f{\Q{s)\ 2 +\M- l P{s)\ 2 )ds+ f\\P{s)\ 2 + \VU{Q{s))\ 2 )ds 
Jo Jo 

The uniform Lipschitz condition Ul implies a linear growth condition on the 
potential force. Hence, there exists a constant C > such that 

\Z(t)\ 2 < \x\ 2 + C f \Z(s)\ 2 ds 
Jo 

By Gronwall's lemma it follows that, 

\Z(t)\ 2 <\x\ 2 e Ch 

for t < h. Hence, for h small enough we obtain the desired bound on the 
second moment of the exact splitting. 
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The difference between f l4"5j) and is, 



Q{h)-Q{h) = 

hM~\l -e^ M ~ lh )P 



r h 

+ / M- 1 [Vt/(Q(s)) - VU(Q(s))](h - s)ds 
Jo 

+ / M- l [- 1 P{s)ds + v / 27/3" 1 rfVF(s)](/i - s) - hM^rj (49) 
Likewise, the difference between (jHj) and ( 14T1) is, 



P{h)-P{h) 



h 



(h-s) 







+ / (e-TM-^ft-s) _ i)vU(Q(s))ds (50) 

From (|49|) and (|50p . it is clear that the leading term of the expectation of these 
differences is 0(h 2 ) and the leading term in the mean-squared expectation 
of the differences is 0(h 3 ^ 2 ). To bound these terms one needs the bounds on 
the second moments of the solutions and the exact splitting provided in (148]) . 
To enable estimation of ( 150]) one needs control of the Hessian of U. The 
assumption of smoothness on U and the uniform Lipschitz condition Ul on 
the potential force provide this control. In particular, since a differentiable 
function is Lipschitz continuous if and only if it has bounded differential, the 
Frobenius norm of the Hessian of U is bounded by the Lipschitz constant of 
the potential force. 

□ 



7.2 Variational Integrators 

Let L : R 2n — >■ R denote the Lagrangian obtained from the Legendre trans- 
form of the Hamiltonian H, and given by: 

L(q,v) = ^v T Mv-U(q). 

A variational integrator is defined by a discrete Lagrangian : M" x M. n x 
R + — > M which is an approximation to the so-called exact discrete Lagrangian 
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which is defined as: 



L$(q , qi ,h)= / L(Q,Q)dt 
Jo 

where Q(t) solves the Euler-Lagrange equations for the Lagrangian L with 
endpoint conditions Q(0) = q and Q(h) = q 1 . 

By passing to the Hamiltonian description, a discrete Lagrangian deter- 
mines a symplectic integrator on R 2n as follows. Given (q ,p ) G M? n , a 
variational integrator defines an update (<fi,Pi) € M. 2n by the following sys- 
tem of equations: 

p = -DtLaiq^q^h), 
p 1 = D 2 L d (q 0} q 1} h). 

Denote this map by h : 1R 2 " -> R 2n , i.e., 

h' (Qo,Po) ^ (Qi,Pi), 

where (qx,Pi) solve (|5T|) . One can show that 8h preserves the canonical 
symplectic form on IR 2ri , and hence, is Lebesgue measure preserving [9]. By 
appropriately constructing Lj, the map Oh can define an approximation to 
the flow of Hamilton's equations for the Hamiltonian H Q. Hyperregularity 
of the discrete Lagrangian means for all h > 

\ det D 12 L d {q , q 1} h)\ > 0, V q , q 1 G Q x Q. (52) 
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